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ABSTRACT 

A comparison of the accuracy of several tech- 
niques recently developed for solving stiff differen- 
tial equations is presented. The techniques examined 
include two general-purpose codes EPISODE and LSODE 
developed for an arbitrary system of ordinary differ- 
ential equations, and three specialized codes CHEMEQ, 
CREK1D and GCKP84 developed specifically to solve 
chemical kinetic rate equations. The accuracy com- 
parisons are made by applying these solution proce- 
dures to two practical combustion kinetics problems. 
Both problems describe adiabatic, homogeneous, gas- 
phase chemical reactions at constant pressure, and 
include all three combustion regimes: induction, heat 

release and equilibration. The comparisons show that 
LSODE is the most efficient code - in the sense that 
it requires the least computational work to attain a 
specified accuracy level - currently available for 
chemical kinetic rate equations. An important finding 
is that an iterative solution of the algebraic enthal- 
py conservation equation for the temperature can be 
more accurate and efficient than computing the tem- 
perature by integrating its time derivative. 

NOMENCLATURE 

ATOL absolute error tolerance for species mole 
numbers 

CPU total computer time required on IBM 370/3033 
computer, sec 

c p 1 constant-pressure molal-specif ic heat of 
species l, J/kmol K 

Erms mean integrated global root-mean-square 
error (Eq. (7)) 


*NRC-NASA Research Associate; on leave from the 
University of Michigan, Dept, of Mechanical 
Engineering and Applied Mechanics, Ann Arbor, 
Michigan 48109. 


EPS for all methods, except EPISODE, local relative 
error tolerance; for EPISODE: local relative 

error tolerance for species with initially non- 
zero mole numbers and for the temperature, and 
local absolute error tolerance for species with 
initially zero mole numbers. 

ERMAX relative error tolerance for Newton-Raphson 
iteration for temperature 

e^ error in ith species mole fraction, (Eq. (4)) 

e rms root-mean-square error in species mole 
fractions and temperature (eq. (6)) 

ey error in temperature (Eq. (5)) 

f 1 net rate of formation of species l (Eq. (1)), 
kmole i / kg mixture sec 

HO initial steplength to be attempted by integra- 
tor, sec 

h 1 molal-specif ic enthalpy of species i, J/kmol 

h 0 mass-specific enthalpy of mixture, J/kq 

NS total number of distinct chemical species in 
mixture 

n 1 mole number of species i, kmole i/kg mixture 

T temperature, K 

T$y standard solution for temperature, K 

t time, sec 

Y-, mole fraction of species l 

Y-, ct standard solution for the mole/f raction of 
’ species i 
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INTRODUCTION 

Many practical chemical reaction flow problems 
require the simultaneous solution of systems of cou- 
pled first-order ordinary differential equations 
(ode's). These ode's describe the time rate of change 
of species concentration and temperature. For species 
l (l = 1,NS) the governing ode can be written as 

dn i 

dt 1 = f 1 ( n k .T); i.k-1. NS 


The motivation for this work is the increasing 
interest in both multidimensional modeling of chemi- 
cally reacting flows and in developing detailed reac- 
tion mechanisms for the combustion of fuels and 
pollutant formation and destruction. Computational 
speed is of primary concern in the former application 
and moderate accuracy is adequate (10-12) . However, 
for developing and validating reaction mechanisms 
accuracy is of critical importance. 

EVALUATION OF TEMPERATURE 


n i (t=0) = given (1) 

T(t=0) = given 

where, n 1 is the mole number of species l (kmole 
l/kg mixture); t is the time (s), T is the tempera- 
ture (K); f } is the net rate of formation of spe- 
cies l (kmole l /kg mixture/sec) due to all forward 
and reverse reactions in which species i partici- 
pates; and NS is the total number of distinct chem- 
ical species in the gas mixture. 

The initial value problem is to solve the system 
of Eqs. (1) for the chemical composition (i.e., n^ 
l » 1, NS) and temperature at the end of a specified 
time interval, given the initial conditions and the 
reaction mechanism. Classical methods such as the 
popular explicit Runge-Kutta method require prohibi- 
tive amounts of computer time to solve large sets of 
chemical kinetic rate equations (_l-4). This is due 
to the extremely small steplengths that classical 
methods have to use to satisfy stability requirements. 
Stable differential equations that impose severe step- 
length limitations on numerical integration routines 
are classified as stiff differential equations (e.g., 
5,6). 


Several solution techniques have been proposed 
and developed for stiff systems of ode's. In Part I 
of this effort and other recent publications ( J.-4 ) , 
the computational work required by several recently 
developed routines in solving chemical kinetic rate 
equations has been examined in detail. In the present 
paper we compare their accuracy. The techniques ex- 
amined in these studies include the general-purpose 
packages EPISODE and LSODE (2-j?), developed as multi- 
purpose stiff differential equation solvers, and the 
specialized methods CHEMEQ (10), CREK1D (21,12) and 
GCKP84 (H,!4), all of which - Rave been developed spe- 
cifically for chemical kinetics applications. These 
solution procedures are summarized in Table I and 
discussed in more detail in Ref. ( 3 ). 

All of the above numerical methods are step-by- 
step methods (i.e., they compute approximations to the 
exact solutions of the ode's at discrete points in 
time). They limit the estimated local truncation er- 
ror (i.e., estimates of the error incurred over one 
time step) to be less than a user-supplied local error 
tolerance. However, the quantity that is of interest 
to the user is the global (i.e., actual) error in- 
curred by a technique in solving the problem. In the 
present paper, which is based on a recent critical 
analysis of the accuracy of the above techniques (15) , 
we provide an estimate of the mean global error and 
examine its variation with the user-supplied local 
error tolerance (EPS) for two practical combustion 
kinetics problems. We also study the computational 
cost (expressed as the computer time required) asso- 
ciated with attaining desired accuracy levels. 


The routines GCKP84 and CREK1D have been devel- 
oped specifically for nomsothermal combustion rate 
equations and therefore include calculation procedures 
for the temperature. For the other techniques, how- 
ever, the temperature had to be solved for along with 
the composition. This was done using one of two dif- 
ferent methods (A and B) described below. 

In the present work, as in Part I, attention is 
restricted to adiabatic, homogeneous, gas-phase chem- 
ical reactions at constant pressure. For such reac- 
tions, the following enthalpy conservation equation 


NS 


L 

i=i 


n^h i = h Q = constant 


( 2 ) 


is an algebraic constraint on the species rate equa- 
tions. In Eq. (2), h^ is the molal-specific en- 
thalpy of species i (J/kmol i) and h 0 is the 
mass-specific enthalpy of the gas mixture ( J/kg) . In 
method A the temperature was calculated from the solu- 
tion for the species mole numbers, n, (i = 1, NS) 
using the initial mixture mass-specific enthalpy, 

Eq. (2) and a Newton-Raphson iteration technique with 
a user-supplied local relative error tolerance, ERMAX. 
The temperature was therefore not an explicit depend- 
ent variable and the integrator tracked only the solu- 
tion for species mole numbers. 

In method B the temperature was evaluated by 
solving its time-derivative obtained by differentiat- 
ing Eq. (2) with respect to time, t 


dT 

dt 


NS 

5 f ' h 


l 


NS 

E 

i=i 


n i c p,i 


(3) 


where Cp -| is the constant-pressure molal 
specif ic-fieat of species i (J/kmol K). In this 
method the temperature was an explicit dependent vari- 
able and the integrator tracked the solutions for both 
the species mole numbers and the temperature. 

TEST PROBLEMS 

The accuracy of the techniques summarized in 
Table I was examined by application to two practical 
combustion kinetics problems. Both problems describe 
adiabatic, homogeneous, constant-pressure, transient, 
batch chemical reaction and include all three combus- 
tion regimes: induction, heat release and 

equilibration. 

Test problem 1 describes the ignition and subse- 
quent combustion of a mixture of 33 percent carbon 
monoxide and 67 percent hydrogen with 100 percent 
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theoretical air at 10 atm and 1000 K initial tempera- 
ture. It involves 12 reactions among 11 species. 

Test problem 2, involving 30 reactions and 15 species, 
describes the ignition and subsequent combustion of a 
stoichiometric hydrogen-air mixture at 2 atm and 
1500 K initial temperature. The reaction mechanisms 
and rate constants for the two problems are given in 
Ref. (3). 

Figures 1 and 2 present the variation of the 
species mole fractions and temperature with time for 
test problems 1 and 2, respectively. Both problems 
were integrated up to time t = 1 ms to obtain near- 
equilibration of all chemical species and temperature. 

COMPUTATIONAL PROCEDURE 

The aim of the present work was to compare the 
accuracy of currently available techniques in solving 
chemical kinetic rate equations. In the absence of 
exact solutions to the test problems, the only way to 
assess the accuracy of a technique is by comparing 
solutions generated with a particular value for the 
local error tolerance (EPS) with those generated with 
a reduced EPS, using either the same technique or a 
different one. To prevent any bias in favor of one 
technique, global errors for each technique were esti- 
mated by comparing its solutions with those generated 
by itself using a reduced EPS. For each technique and 
test problem, the solution used as a basis of compar- 
ison was the most accurate generated and is referred 
to as its standard solution. For example, for CREKlD 
standard solutions were generated with CREKlD and EPS 
= 10~% These standard solutions were compared with 
runs using CREKlD and EPS = 10” 2 , 10% and lO -4 
to assess the global errors incurred by the latter. 


mole fractions and temperature were generated. De- 
tailed plots of e-j(t) and ej(t) are presented 
in Ref. (lj[) . 

Standard Solutions . For consistency, standard 
solutions were generated with the same value of EPS 
(=10 -5 ) for all methods, except EPISODE. In con- 
trast to the other routines for which EPS is the local 
relative error tolerance, it is a mixed error toler- 
ance for EPISODE - relative for species with initially 
nonzero mole numbers (i.e., reactants) and for the 
temperature; and absolute for species with initially 
zero mole numbers (i.e., all intermediate species and 
products) (15) . EPISODE is therefore inferior to the 
other codes - i.e., for a given value of EPS, it pro- 
duces less accurate solutions (3,lj[). With EPISODE, 
we have therefore used the smallest value of EPS 
(=1CT 6 ) that the single-precision version allowed. 

All the techniques examined here required the 
specification of several input parameters, in addition 
to the local error tolerance, EPS. In generating 
standard solutions for each technique, the values used 
for the input parameters were those that resulted in 
the most accurate solution - see Ref. (15) for de- 
tails. For 6CKP84, however, since detaTTs of the in- 
tegration technique are not yet known, default values 
were used for all parameters except the initial step- 
length (HO) to be attempted by the ode solver. In our 
previous work (1-3) with GCKP84, a default value of 
HO = 10’ 6 s hacT Eeen used. However, Bittker and 
Scullin (14) have since then set the default value for 
HO at 5xTO“ 8 sec. We have, nevertheless, used HO 
= 10* 6 sec because it produced more accurate results. 

RESULTS AND DISCUSSION 


A typical computational run consisted of initial- 
izing the time (t, set equal to zero), the species 
mole numbers and the temperature. The integrator was 
called with values for the necessary input parameters 
including the local error tolerance (EPS) and the time 
at which the integration was to be terminated (= 1 ms 
for both problems). The values used for the other 
input parameters, obtained in a previous study (_3) by 
a tri al-and-error procedure, resulted in the least CPU 
time for each technique and value of EPS. After each 
step successfully executed by the code, values for the 
time reached by the integrator and the species mole 
fractions and temperature computed at that time were 
saved. The values saved for the time served as input 
data for the output stations at which the standard 
solution was to be generated. The global errors in 
the solutions for the species mole fractions and tem- 
perature were estimated by comparisons with the stan- 
dard solution values for these quantities as follows. 


e^t) = y 


Y, (t) 


i,ST 


rti 


- 1 


e T (t) = I 


T(t) 


ST 


(tT 


(4) 

(5) 


In Eqs. (4) and (5) e.,(t) and e-j-(t) are, re- 

spectively, the global errors at time t, in the mole 
fraction, Vi(t), of species i and the temperature, 

T ( t ) , Y-, 5 t( t) and T^j ( t ) are, respectively, the 
standard* solution values for the mole fraction of spe- 
cies i and the temperature at time, t. The global 
errors in species mole fractions and temperature and 
the corresponding times were saved for later analysis. 
In this way time histories of the errors in species 


The procedure outlined in the section Computa- 
tional Procedure was used to examine the global errors 
incurred by all techniques in solving the two problems 
described in the section Test Problems. All results 
presented herein were obtained on the NASA Lewis 
Research Center's IBM 370/3033 computer using single- 
precision accuracy, except GCKP84 which was in 
double-precision. 

With LSODE, EPISODE and CHEMEQ, both temperature 
methods A and B were attempted. The following naming 
convention was used. Techniques using temperature 
method A were given the suffix A (LSODE-A, EPISODE-A 
and CHEMEQ-A) and those using temperature method B 
were given the suffix B (LSODE-B, EPISODEB and 
CHEMEQ-B). For consistency between the two methods, 
ERMAX (the maximum local relative error allowed in the 
temperature in method A) was set equal to EPS which is 
approximately the maximum local relative error allowed 
in the temperature in method B. 

To facilitate accuracy comparisons, at each value 
of the time where global errors had been evaluated 
(and saved), the root-mean-square (rms) error, e rms (t) 
was computed using 


e 


Jt) 



e?(t) + e|(t) 
NS + 1 


( 6 ) 


To prevent the possibility of requiring accuracy 
from species with immeasurably small concentrations, 
the summation in Eq. (6) does not include species 
whose standard solution mole fractions were less than 
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0.1 ppm - for such species, e-,(t) was set equal to 
zero. 

The maximum percent root mean square error in- 
curred by all techniques are given in Tables II and 
III, respectively, for test problems 1 and 2. Also 
given in these tables are the values used for the in- 
put parameters required by each method. These values, 
taken from Radhakrishnan resulted in the least 
CPU time (CPU in Tables II and III) to solve the prob- 
lem with the given value of EPS. In these tables, HO 
is the user-supplied value for the initial steplength 
to be attempted by the integrator, and ATOL is the 
absolute error tolerance for all species mole numbers 
- a value of zero was used for the absolute error 
tolerance for the temperature (required by LSODE-B). 


to ensure accuracy can result in excessively large CPU 
times. For example, for test problem 2 the run using 
LSODE-B with EPS = 10“ 5 and ATOL = 10“ n required 
about 3.4 sec CPU time; in contrast, the run with 
ATOL = 1(T 15 required almost 20 sec although the 
solution was not significantly different. In addi- 
tion, as discussed in Part I (4), the solution can be 
adversely affected by a poor clToice for the output 
stations at which the solution is desired. 

Tables II and III show also that the use of tem- 
perature method A does not result in significantly 
larger errors than those incurred by temperature meth- 
od B. On the contrary, method A can be more accurate 
than method B - this difference in accuracy is most 
marked for CHEMEQ. 


For test problem 1, the run with GCKP84 and EPS 
= 10“ 2 displayed serious instability and so was ter- 
minated. For values of EPS > 5x10 , EPISODE-A and 

-B predicted little or no change from the initial com- 
position or temperature after an elapsed time of 1 ms. 
Similar remarks apply to test problem 2 and the runs 
with EPISODE-A and EPS > 5xl0 -4 and those with 
EPISODE-B and EPS >_ 5xl(5“ 2 . Although the runs with 
EPISODE-B and EPS = 10“ 3 and 5xl0“ 4 were successfully 
completed in that correct solutions were returned at 
t = 1 ms, they were significantly inaccurate during 
heat release. For example, the run with EPS = 5xl0“ 4 
predicted little or no change from the initial condi- 
tions until t = 40 ps when heat release was predicted 
to start. In contrast, for the standard solution heat 
release is almost over by this time (Fig. 2). 

Tables II and III show large variations in the 
maximum root mean square error incurred by the dif- 
ferent techniques. GCKP84 and EPISODE experienced the 
greatest difficulty in tracking their standard solu- 
tions and maximum root mean square errors of over 100 
percent were obtained with these codes. In contrast, 
the maximum errors incurred by CHEMEQ, CREKlD and 
LSODE were significantly smaller. Comparisons of the 
runs with the largest values of EPS show that LSODE 
was the most accurate code for test problem 1, and 
CREKlD for test problem 2. 

The accuracy obtained with LSODE was strongly 
influenced by the value selected for ATOL (the abso- 
lute error tolerance for the species mole numbers). 

For problem 2 LSODE-A showed little sensitivity to 
changes in EPS, because all runs used the same ATOL 
(Table III). In contrast, the runs with LSODE-B 
showed decreased errors with reductions in EPS because 
of the use of smaller values of ATOL. Note the sig- 
nificant reduction in the maximum root mean square 
error obtained by decreasing ATOL from 10“° to 10 - * 2 
for LSODE-B and EPS = 10“ 4 (Table III). This de- 
crease in the maximum error with a reduction in ATOL 
is further illustrated for problem 2 by the results 
presented in Table IV for LSODE-B and EPS = 10“^. 

These results show that great care must be exercised 
in specifying values for ATOL. A poor guess for ATOL 
can result in significantly inaccurate solutions. The 
CPU time required by LSODE was also affected by ATOL. 
In general, a trial-and-error procedure was necessary 
to obtain the optimum value for ATOL - defined as that 
value which results in the least CPU time while satis- 
fying prescribed accuracy requirements. This trial- 
and-error search - which can be time consuming, 
especially for large systems of ode's - for the opti- 
mum value for ATOL is the main difficulty associated 
with using LSODE for solving chemical kinetic rate 
equations. The use of extremely low values for ATOL 


To provide a more comprehensive measure of the 
accuracy than the maximum root mean square error, we 
have adopted the following procedure. For each run a 
mean integrated root mean square error, E rms , was 
calculated as follows: 


E 


rms t 


1 

end 



e me (t) dt 
rms 


( 7 ) 


where t enc j is the end of the prescribed time 
interval . 

Equation (7) provides a single quantity that is 
a measure of the average error incurred in solving the 
complete problem. The integral was evaluated numeri- 
cally using Simpson's rule modified for unequal step 
sizes. 


Figures 3 and 4 present the variation of 
with the user-supplied local error tolerance, EPS. In 
addition, Table IV gives values of E rm s incurred by 
the different runs with LSODE-B and EPS = 10“ 5 for test 
problem 2. These results illustrate the increasing 
accuracy obtained with reductions in ATOL and the sig- 
nificant errors that can be incurred by a poor guess 
for ATOL. Figures 3 and 4 show that all methods used 
in the present study are tolerance-effective (i.e., 
decreasing EPS results in reduced error). For both 
problems temperature method A is as accurate as meth- 
od B - in many cases, it is significantly more accu- 
rate. For problem 2 and EPS = 10“ 4 , LSODE-B is more 
accurate than LSODE-A because of the smaller ATOL 
used. 


Figures 3 and 4 illustrate the significant dis- 
crepancies between the values specified for EPS and 
the errors actually obtained. For example, for prob- 
lem 1 a value of EPS = 10“ 2 (1 percent) has resulted 
in an average error of almost 50 percent error for 
CHEMEQ-B. These plots show that for a given value of 
EPS, LSODE is the most accurate code currently avail- 
able for solving chemical kinetic rate equations. 
However, GCKP84, CREKlD and CHEMEQ-A compare very 
favorably with LSODE for small values of EPS. EPISODE 
is significantly less accurate than the other codes 
because the error control used in it is inappropriate 
for chemical kinetics applications (_3). 

Figures 5 and 6 present the variation of the com- 
putational work (expressed as the CPU time in seconds 
required on the NASA Lewis Research Center's IBM 370/ 
3033 computer) with the mean integrated error for test 
problems 1 and 2, respectively. For both problems the 
CPU times required by temperature method A are less 
than, or compare favorably with, those required by 
method B to attain the same accuracy levels. This 
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difference is most pronounced for CHEMEQ and EPISODE. 
Note that for EPISODE— B the computational work in- 
creases with increasing error. 

Figures 5 and 6 illustrate the large differences 
in the computational work required by the different 
techniques to attain comparable accuracy levels. For 
problem 1 and a one-half percent mean global error, 
the CPU time varies from about 0.37 sec for LSODE-A 
to over 40 sec for CHEMEQ-A. For problem 2 this dif- 
ference is even greater. For both test problems LSODE 
is the most efficient code in the sense that it re- 
quires the least CPU time to attain a specified accu- 
racy level. 

For test problem 2 EPISODE-A compares very favor- 
ably with LSODE (Fig. 6). However, the solution gen- 
erated by EPISODE was found to strongly depend on the 
value selected by the user for the initial steplenqth 
(HO) to be attempted by the integrator. A poor guess 
for HO can result in inaccurate and unstable solu- 
tions (1-4). can also result in excessive CPU 
times. For example, the run using EPISODE-A with EPS 
= 10~ 4 and HO = 10“° sec required almost 129 sec 
for problem 2; in contrast, the run with HO = 10” 7 
sec required only about 0.6 sec. 

To attain the same accuracy level, CREK1D 
requires more CPU time than LSODE for both test prob- 
lems. CREK1D is, however, attractive for the follow- 
ing reason. It is intended primarily for performing 
multidimensional calculations of chemically reacting 
flows by coupling it with a hydrodynamic equation 
solver. Currently available hydrodynamic codes are 
at best accurate to within a few percent, so genera- 
tion of highly accurate chemical kinetic solutions is 
wasteful (_10) . For problem 1 CREKlD requires approx- 
imately 0.3 sec to produce a solution with a mean 
global error of about 5 percent, which is sufficiently 
accurate for multipoint calculations of reacting 
flows. This CPU time compares favorably with the 

0.37 sec required by LSODE-A. 

CONCLUSIONS 

A comparison of the accuracy of several recently 
developed numerical techniques (GCKP84, CREKlD, LSODE, 
EPISODE, and CHEMEQ) in solving chemical kinetic rate 
equations has been made. This study has shown that 
LSODE is the most efficient code - in the sense that 
it requires the least CPU time to attain a specific 
accuracy level - currently available for solving chem- 
ical kinetic rate equations. The major difficulty 
associated with its use is the tri al-and-error proce- 
dure necessary to obtain optimum values for the abso- 
lute error tolerances for the variables. A poor guess 
for the absolute error tolerance can result in exces- 
sive CPU times or in seriously inaccurate solutions. 
Wnen accuracy is not of primary concern, as in multi- 
dimensional chemically reacting flow calculations, 
CREKlD is an attractive alternative to LSODE. 

An important conclusion is that the use of an 
algebraic enthalpy conservation equation for calculat- 
ing the temperature can be more accurate and efficient 
than evaluating the temperature by integrating its 
time-derivative. 
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TABLE I - SUMMARY OF SOLUTION PROCEDURES EXAMINED 


Method 

Description 

GCKP84 

Details not yet available 

CREK1D 

Variable-step, predictor-corrector method based on an exponentially 
fitted trapezoidal rule, includes filtering of ill-posed initial 
conditions and automatic selection of functional iteration or Newton 
iteration 

LSODE 

Variable-step, variable-order backward differentiation method with a 

EPISODE 

generalized Newton iteration 3 

CHEMEQ 

Variable-step, second-order predictor-corrector method with an 
asymptotic integration formula for stiff equations 


3 Other options are included in these packages but this option was found to be 
the fastest for both problems, and is the only one considered in this work. 


TABLE II - MAXIMUM PERCENT RMS ERRORS AND MINIMUM 
CPU TIMES REQUIRED FOR TEST PROBLEM 1 


Method 

EPS 

ATOL 

HO, 

sec 

Max rms 

error, 

percent 

CPU, 

sec 

GCKP84 

a 9 


bio -6 

467 9 

(a) 


5xl(H 


bio-6 

467 9 

86 


itH 


biQ-6 

258 6 

.98 


10” 4 


biQ- 6 

115 0 

1.13 

CREK1D 

10-2 



30 79 

27 

lO" 3 



7 218 

33 


10- 4 



2 311 

67 

LSODE-A 

10-2 

10~u 


10.52 

.37 


lO" 3 

10-3} 


8 596 

.46 


10- 4 

ltr 11 


4.786 

.63 

LSODE-B 

10-2 

lO-}? 


10.57 

40 


10-3 

10-11 


8 167 

.44 


10~ 4 

10-11 


6 366 

59 

EPISODE-A 

5xio- 6 


lO"? 

825.7 

.034 

EPISODE-B 

5xl0~ 6 


r-. 

1 

O 

825.9 

.035 

CHEMEQ-A 

10-2 



34.47 

15.1 


10“^ 



6.580 

28.4 


10 -4 



.868 

39 8 

CHEMEQ-B 

10-2 



57.61 

15.5 


1 o -4 



19.48 

22 7 


10 -4 


— 

2.219 

36 5 


terminated at t = 3 7x10“^ sec because of problems 

with serious instability. 

^Default value 





TABLE III. - MAXIMUM PERCENT RMS ERRORS AND MINIMUM 
CPU TIMES REQUIRED FOR TEST PROBLEM 2 


Method 

EPS 

ATOL 

HO, 

sec 



GCKP84 

10-| 


a io-6 

227.2 

1.85 


icr 3 


aio-6 

98.65 

1.91 


10" 4 

— 

a 10 -6 

21.76 

2.44 

CREK1D 

10-2 



3 959 

1 18 

icr 3 



2 148 

1 07 


icr 4 

— 

— 

0.612 

2.32 

LSODE-A 

10-2 

icr 5 


29.38 

.54 


icr 3 

icr 8 


27.37 

.78 


icr 4 

icr 8 


29.44 

.94 

LSODE-B 

10-2 

10- 3 


30 53 

.52 


10" 3 

IO" 9 


3.484 

.94 


icr 4 

10-8 


1 528xl0 3 

1.49 


io- 4 

IO" 12 


8.641 

1.46 

EPISODE-A 

5xicr 4 


10-7 

290 4 

065 


icr 4 


10-7 

121 8 

.59 


icr 5 


IO -6 

54 08 

.78 

EPISODE-B 

5xicr 4 


10-10 

280.1 

2.38 


10- 4 


IO" 11 

131 6 

91 


icr 5 


icr 9 ’ 

78.31 

.88 

CHEMEQ-A 

10~2 



26.97 

77 7 

10r3 



2.133 

mm 


IO" 4 

— 


.207 

76.1 

CHEMEQ-B 

10-2 





23 80 

36.3 


IO -3 





4.631 

43.0 


IO" 4 


— 

2.341 

87.6 


a Default value 


TABLE IV. - EFFECTS OF ABSOLUTE ERROR 
TOLERANCE FOR SPECIES MOLE NUM- 
BERS (ATOL) ON MAXIMUM PERCENT 
ROOT MEAN SQUARE ERROR, 

MEAN INTEGRATED ROOT 
MEAN SQUARE (E rm$ ) 

AND CPU TIME 

[All results obtained with LSODE-B 
and EPS=10” 3 for probl em 2.] 


ATOL 

Max rms 
error, 
(percent ) 

^rms 

CPU, 

(sec) 

10~ 8 

1 112xl0 4 

2.792X10" 1 

4 17 

10 ”?n 

2 022xl0 3 

4.011x10-2 

3 24 

10-!° 

20.25 

9.421xl0~ 4 

2.52 












10~ 6 10- 5 10' 4 1 CT 3 

Time, t, s 


Figure 1. - Variation of temperature and species 
mole fractions with time for test problem 1. 
Solution generated using LSODE-B with EPS = 
IQ -5 . 




10' 6 10" 5 10' 4 10" 3 


Time, t, s 

Figure 2, - Variation of temperature and species 
mole fractions with time for test problem 2. 
Solution generated using LSODE-B with EPS = 
10‘ 5 . 






Figure 3. - Variation of the mean integrated global error (E rms ) 
with the local error tolerance (EPS) for test problem 1. 



MEAN INTEGRATED ERROR. 








-2 


10 
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D ERROR, E rms 


with the mean integrated global error (E rms ) 
370/3033 computer. For symbol key, see 
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